Energy transport in disordered classical spin chains 
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We present a numerical study of tlie diffusion of energy at higii temperature in strongly disordered 
chains of interacting classical spins evolving deterministically. We find that quenched randomness 
strongly suppresses transport, with the diffusion constant becoming reduced by several orders of 
magnitude upon the introduction of moderate disorder. We have also looked for but not found signs 
of a classical many-body localization transition at any nonzero strength of the spin-spin interactions. 
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I. INTRODUCTION 

Chemical, structural and other imperfections are un- 
avoidable in a typical crystalline solid. Often, reducing 
their concentration reveals interesting intrinsic proper- 
ties of the ideal "clean" material. It has long been recog- 
nized, however, that sufficiently strong disorder can lead 
to a host of phenomena essentially distinct from what is 
observed in clean materials. Anderson localization fl| is 
perhaps the best known example of such a phenomenon, 
whereby a quantum particle becomes coherently trapped 
by the impurity potential and no longer participates in 
transport. Essentially a wave phenomenon, Anderson lo- 
calization can also occur for classical linear waves, e.g. 
photons or phonons. 

It has often been assumed that true localization 
(strictly zero diffusivity) of interacting particles can only 
happen at absolute zero temperature, even though An- 
derson's original paper on localization discusses the pos- 
sibility of localization persisting at nonzero temperature 
Recently this question has been examined carefully 
by Basko and collaborators 0, who performed a stabil- 
ity analysis of an Anderson insulator against weak in- 
terparticle interactions at low but nonzero temperature. 
Their central conclusion is that an isolated system of 
strongly disordered but weakly interacting quantum par- 
ticles should exhibit a transition into an insulating phase 
with strictly zero diffusion at some low but nonzero ex- 
citation energy per particle (or temperature) . Motivated 
by this work, two of us considered quantum lattice 
models that would be expected to exhibit such a dy- 
namical many-body localization transition as one varies 
the interactions or disorder strength, even at arbitrarily 
high temperature. We attempted to detect this tran- 
sition by studying exact many-body spectral statistics 
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of small lattices. Initial results are encouraging fs'], al- 
though they fall short of making a strong case for the 
existence of a quantum many-body localization transi- 
tion, due to strong finite-size effects. 

Setting aside the question of the existence of such a 
transition (i.e., assuming it does exist), one might won- 
der about its nature, e. g., the universality class. On 
the one hand, the theoretical analysis of Basko, et al. 
relies entirely on quantum many-body perturbation the- 
ory. Rather generally, however, one expects macroscopic 
equilibrium and low-frequency dynamic properties of in- 
teracting quantum systems at nonzero temperature to be 
describable in terms of effective classical models. This ex- 
pectation is certainly borne out in a variety of symmetry- 
breaking phase transitions with a diverging correlation 
length, such as, e. g., a finite temperature Neel ordering 
of spin-1/2 moments. One can begin to understand the 
microscopic mechanism behind such a many-body "cor- 
respondence principle" as a consequence of an effective 
coarse-graining, whereby the relevant degrees of freedom 
are correlated spins moving together in patches that grow 
in size as the phase transition is approached and there- 
fore become "heavy" and progressively more classical. 
Further extension of these ideas to general, non-critical, 
dynamical response is more involved: roughly speaking, 
it requires that the typical many-body level spacing in 
each patch be much smaller than the typical matrix el- 
ement of interactions with other patches. If this is true 
(as it is in most models at finite temperature, though 
not necessarily in the insulating phase analysed by Basko 
and collaborators) one replaces microscopic quantum de- 
grees of freedom with macroscopic classical ones, which 
typically obey "hydrodynamic" equations of motion at 
low frequencies [J|. Since it is expected that the many- 
body localization transition is accompanied by a diverg- 
ing correlation length (akin to the Anderson transition) 
one might expect some sort of classical description to 
emerge en-route from the localized phase to the diffusive 
phase. It was this thinking that initially motivated us to 
consider the possibility of classical many-body localiza- 
tion. 
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The process by which cohective classical (hydro-) dy- 
namics emerges from a microscopic quantum description 
is subtle and may or may not be relevant to the many- 
body localization discussed above. A somewhat less sub- 
tle, but apparently largely unexplored related question, 
is whether nonlinear, interacting and disordered classi- 
cal many-body systems are capable of localization at 
nonzero temperature. To be precise, a many-body classi- 
cal dynamical system with a local Hamiltonian (including 
static randomness) should show hydrodynamic behavior, 
e.g. energy diffusion, provided the local degrees of free- 
dom are nonlinear and interacting, and the disorder is not 
too strong. In this regime, the isolated system can func- 
tion as its own heat bath and relax to thermal equilib- 
rium. Diffusive energy transport must stop if the interac- 
tions between the local degrees of freedom are turned off. 
How is this limit approached? Can there be a classical 
many-body localization transition where the energy dif- 
fusivity vanishes while the interactions remain nonzero? 
These are the basic questions we set out to investigate in 
this paper. 

Our preliminary conclusion is that classical many-body 
systems with quenched randomness and nonzero nonlin- 
ear interactions do generically equilibrate, so there is no 
generic classical many-body localized phase. Our picture 
of why this is true is that generically a nonzero frac- 
tion nonlinearly interacting classical degrees of freedom 
are chaotic and thus generate a broad-band continuous 
spectrum of noise. This allows them to couple to and 
exchange energy with any other nearby degrees of free- 
dom, thus functioning as a local heat bath. Random 
classical many-body systems generically have a nonzero 
density of such locally-chaotic "clusters", and thus the 
transport of energy between them is over a finite dis- 
tance and can not be strictly zero, resulting in a nonzero 
(although perhaps exponentially small) thermal conduc- 
tivity. Quantum systems, on the other hand, can not 
have a finite cluster with a truly continuous density of 
states: the spectrum of a finite cluster is always discrete. 
Thus the mechanism that we propose forbids a generic 
classical many-body localized phase, yet it does not ap- 
pear to apply to the quantum case. The proposed exis- 
tence of the many-body insulator in quantum problems 
is then a remarkable manifestation of quantum physics in 
the macroscopic dynamics of highly-excited matter. In 
this paper we shall primarily focus on macroscopic low 
frequency behavior, postponing detailed analysis of local 
structure of noise and its relation to transport. Our con- 
clussions are broadly consistent with findings of Dhar and 
Lebowitz[6] although given the rather major differences 
in models, methods and, most importantly, the extend 
to which strongly localized regime is probed we refrain 
from making direct comparisons. 

We study energy transport in a simple model of local 
many-body Hamiltonian dynamics that has both strong 
static disorder and interactions: classical Heisenberg spin 
chains with quenched random fields. For simplicity, we 
consider the limit of infinite temperature, defined by av- 



eraging over all initial conditions with equal weights. Our 
systems conserve the total energy and should exhibit en- 
ergy diffusion; they have no other conservation laws. The 
energy diffusion coefficient, D, can be deduced from the 
autocorrelations of the energy current (as explained be- 
low) and is shown in Fig. [T] as a function of the strength 
of the spin-spin interactions, J. The mean-square ran- 
dom field is A-^, and as we vary J we keep 2J^ + = 1, 
as explained below. The limit J is where the inter- 
actions vanish, so there is (trivially) no energy transport. 
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FIG. 1: Disorder-averaged energy diffusion constant D as a 
function of the spin-spin interaction J. The line has slope 8 
on this log-log plot. 

As the interaction J is decreased, the thermal diffu- 
sivity D decreases very strongly; we have been able to 
follow this decrease in D for about 5 orders of magnitude 
before the systems' dynamics become too slow for our nu- 
merical studies. For most of this range, we can roughly 
fit D{J) with a power-law, D J'^, with a rather large 
exponent, 7 = 8, as illustrated in Fig. 1. This large 
exponent suggests that the asymptotic behavior at small 
interaction J may be some sort of exponential, rather 
than power-law, behavior, consistent with the possibility 
that the transport is actually essentially nonperturbative 
in J. In principle, it is also possible to fit these data to 
a form with a nonzero critical Jc, so that D{J < Jc) — 
— such fits prove inconclusive as they produce estimates 
of Jc considerably smaller than the values of J where we 
can measure a nonzero D. Since we are not aware of any 
solid theory for the behavior of D{J), these attempts at 
fitting the data are at best suggestive. The large range of 
variation of the macroscopic diffusion constant D across 
a rather modest range of J is the most clearly remarkable 
and robust finding that we wish to present in this paper. 

Our model and general methods employed will be pre- 
sented and discussed in the next Section. Much of what 
we present is based on the analysis of energy current fiuc- 
tuations in isolated rings. For various reasons we have 
found it beneficial to focus on these rather than fiuctua- 
tions of the energy density or on current carrying states 
in open systems (we have spot-checked for quantitative 
agreement among these three methods). In Section HI 
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we present our results for macroscopic transport starting 
from short time behavior that is relatively easy to un- 
derstand and working up to long time, DC behavior that 
is both difficult to compute and as of yet poorly under- 
stood. One particularly interesting observation we make 
here is that of a subdiffusive behavior at long times, ap- 
parently distinct from the much discussed mode-coupling 
behavior well representative of linear diffusion in the 
presence of disorder. We discuss some afterthoughts and 
open problems in the Summary, with some important 
additional details relegated to the Appendices (such as 
quantitative explorations of finite-size effects, roundoff, 
many-body chaos and self-averaging). 

II. MODEL, TRAJECTORIES AND 
TRANSPORT 

The classical motion of N interacting particles is usu- 
ally defined by a system of coupled differential equations 
of motion. The "particles" we study here are classical 
Heisenberg spins - three-component unit-length vectors. 
Si, placed at each site i of a one-dimensional lattice. 
With a standard angular momentum Poisson bracket and 
a Hamiltonian, H , the equations of motion are 

§^-H,xS,, (1) 

where Hj = dH/dSj is the total instantaneous field act- 
ing on spin Sj. The Hamiltonians we consider are all of 
the form 

if = ^(h, .8, + JS,.S,+i), (2) 

j 

with uniform pairwise interaction J between nearest- 
neighbor spins, and quenched random magnetic fields, 
hj. For almost all of the results in this paper, we choose 
the random fields to be hj = hjfij, where the hj are 
independent Gaussian random numbers with mean zero 
and variance A^, while the iij are independent randomly- 
oriented unit vectors, uniformly distributed in orienta- 
tion. Because of the random fields, total spin is not 
conserved and we can focus on energy diffusion as the 
only measure of transport in this system. For J — 
and A > any initial distribution of energy is localized, 
as the spins simply precess indefinitely about their local 
random fields, so the diffusivity is -D = 0. In the opposite 
limit, where A = and J > 0, there is diffusive trans- 
port with D ^ J (with nonlinear corrections due to the 
coupling between energy and spin diffusion [1, Q ) • We 
are interested in the behavior of D as one moves between 
these two limits, especially as one approaches J = with 
A > 0. 

Given initial spin orientations, it is in principle 
straightforward to integrate the equations of motion nu- 
merically, thus producing an approximate many-body 
trajectory. Correlation functions can then be computed 



and averaged over a such trajectories and over realiza- 
tions of the quenched random fields. The transport co- 
efficients can thereby be estimated via the fluctuation- 
dissipation relations. 



A. The Model 

Before we embark on this program, however, we start 
by making a change to the model's dynamics ([1]), but 
not to its Hamiltonian in order to facilitate the nu- 
merical investigation of the long-time regime of interest 
to us, where the diffusion is very slow. In order to get 
to long times with as little computer time as possible, 
we want our basic time step to be as long as is possible. 
What we are interested in is not necessarily the precise 
behavior of any specific model, but the behavior of the 
energy transport in a convenient model of the type 
Since we are studying energy transport, it is absolutely 
essential that the numerical procedure we use does con- 
serve total energy (to numerical precision) and that the 
interactions and constraints remain local. Thus we mod- 
ify the model's dynamics to allow a large time step while 
still strictly conserving total energy. 

We change the equations of motion ([1]) of our model so 
that the even- and odd-numbered spins take turns pre- 
cessing, instead of precessing simultaneously. We will 
usually have periodic boundary conditions, so we thus 
restrict ourselves to even length (thus bipartite) chains. 
We use our basic numerical time step as the unit of time 
(and the lattice spacing as the unit of length) . During one 
time step, first the odd-numbered spins are held station- 
ary, while the even-numbered spins precess about their 
instantaneous local fields, 

Hrit) ^ hr + JSr-l{t) + JSr+l{t) , (3) 

by the amount they should in one unit of time accord- 
ing to (dl). Note that since the odd spins are station- 
ary, these local fields on the even sites are not changing 
while the even spins precess, so that this precession can 
be simply and exactly calculated, and the total energy 
is not changed by this precession. Then the even spins 
are stopped and held stationary in their new orienta- 
tions while the odd spins "take their turn" precessing, to 
complete a full time step. Although this change in the 
model's dynamics from a continuous-time evolution to a 
discrete-time map is substantial, we do not expect it to 
affect the qualitative long-time, low-frequency behavior 
of the model that is our focus in this paper. In partic- 
ular, we clearly observe correct diffusive decay of local 
correlations for weak disorder and essentially indefinite 
precession of spins at very strong disorder. 

We have decided to use parameters so that the mean- 
square angle of precession of a spin during one time step is 
one radian (at infinite temperature), which seems about 
as large as one can make the time step and still be roughly 
approximating continuous spin precession. This choice 
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dictates that the parameters satisfy 



2J^ 



= 1 



(4) 



We will generally describe a degree of interaction by quot- 
ing the J; the strength A of the random field varies with 
J as dictated by 



B. Observables 

The basic observable of interest, the instantaneous en- 
ergy ei{t) at site i is 



e^it) 



h, -8,(0 + ^(8, 



-i(<)-S,(i) + S,+i(t).S,(i)) . (5) 



Note that with this definition, the interaction energy cor- 
responding to a given bond is split equally between the 
two adjacent sites. When updating the spin at site i, 
only the energies of the three adjacent sites, and ei±i, 
change, due to the change in the interaction energies in- 
volving spin i. This rather simple pattern of rearrange- 
ment of energy allows for an unambiguous definition of 
the energy current at site i during the time step from 
time t to t + 1. If site i is even, so it precesses first, then 
the current is 

j,(t) = J[S,{t + 1) - S,{t)] ■ [S,+i(t) - S,_i(<)] , (6) 

while for i odd, 

Mt) - ./[S,(t + l)-S,(t)].[S,+i(i+l)-S,_i(< + l)] . (7) 

We are working at infinite temperature, or alterna- 
tively at /3 = {ksT)^^ = 0. The conventionally-defined 
thermal conductivity vanishes for /3 ^ [1]. Instead, 
here we define the DC thermal conductivity n so that 
the average energy current obeys 



(8) 



in linear response to a spatially- and temporally-uniform 
small gradient in /? = \/{kBT). The Kubo relation then 
relates this thermal conductivity at /3 = to the correla- 
tion function of the energy current via 



where 



c{t) ^Y.^{jmm))] 



(9) 



(10) 



is the autocorrelation function of the total current, where 
the square brackets, [...], denote a full average over in- 
stances of the quenched randomness ( "samples" ) and the 
angular brackets, (...}, denote an average over initial con- 
ditions in a given sample and time average within a given 
run. For our model ^ the average energy per site obeys 



dm 

dp 



J2 



(11) 



at /3 = 0, and the energy diffusivity D is then obtained 
from the relation 



D 



dm 

d(3 



(12) 



In a numerical study, if a quantity (such as k) is non- 
negative definite, then it is helpful to measure it if possi- 
ble as the square of a real measurable quantity. We use 
this approach here, noting that 



(13) 



For a particular instance of the random fields in a chain 
of even length L with periodic boundary conditions and 
a particular initial condition / run for time t, we thus 
define the resulting estimate of k as 



(14) 



If these estimates are then averaged over samples and 
over initial conditions for a given L and this results 
in the estimate KL{t) = [(K/(i))] . These estimates ^^{1) 
must then converge to the correct DC thermal conduc- 
tivity K, in the limits i, t ^ oo. 



C. Finite-size and finite-time eff'ects 

In a sample of length L, we expect finite-size effects to 
become substantial on time scales 



t> II ^ CdL^ /Deff, 



(15) 



where Dctf is the effective diffusion constant at those time 
and length scales, and we find Co = 10 (remarkably Eq. 
15 remains valid more or less with the same value of Cjj 
across the entire range of parameters - see Appendix [X| . 
With periodic boundary conditions (which is the case 
in our simulations) this means that KL(t) saturates for 
t > t-L to a, value different from (and usually above) its 
true DC value in the infinite L limit, while with open 
boundary conditions (no energy transport past the ends 
of the chain) the infinite-time limit of KL{t) is instead 
identically zero for any finite L. We simply avoid this 
purely hydrodynamic finite-size effect by using chains of 
large enough length L, which is relatively easy, especially 
in the strongly-disordered regime of interest, where D^ff 
is quite small. 

For the smallest values of J that we have studied, the 
system is essentially a thermal insulator, and the D^jf 
is so small that finite-size effects are just not visible at 
accessible times even for small values of L, such as L = 
10. Instead, given the way we are estimating k, a finite- 
time effect, due to the sharp "cutoffs" in time at time 
zero and t in dominates the estimates K^it) ^ J"^ /t 
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in this small-J regime. To explain this better we can 
rewrite the definition of k as 

t t t/2 

'^^ W = II E E - = ;^ E ^L{ta.\ (16) 

t=1t' = 1 tav = l 

where we have assumed an even t (there is an additional 
term otherwise) and k\{t) = C{t'). Localization, 
i.e. zero DC conductivity, implies a rapidly vanishing k* 
as well as k at long times. The latter however acquires 
a tail, K ~ 1/i whose amplitude is set by the short-time 
values of k*. 

For the intermediate values of J that are of the most 
interest to us in this paper, there is also another, stronger 
finite-time effect due to an apparently power-law "long- 
time tail" in the current autocorrelation function, C(r), 
as we discuss in detail below. Importantly, at long times 
this intrinsic finite-time effect dominates the extrinsic, 
cutoff- induced, 1/t effect discussed above, so KL{t) re- 
mains a useful quantity to study in this regime. 

III. RESULTS: MACROSCOPIC DIFFUSION 
A. Current autocorrelations 

Since the total current is not dynamically stationary, 
its autocorrelation function, C{t), should decay in time. 
In a strongly disordered dynamical system we expect the 
DC conductivity, which is the sum over all times of this 
autocorrelation function, to be very small due to strong 
cancelations between different time domains (i.e., C'(t) 
changes sign with varying t). The basic challenge of 
computing the DC thermal conductivity k boils down 
to computing (and understanding) this cancelation. 

The autocorrelation function C{t) has three notable 
regimes as we vary J and t. First, C{t) is positive and 
of order at times less than or of order one, as illus- 
trated in Fig. [21 It quickly becomes negative at larger 
times. For small J it is negative and of order in mag- 
nitude for times of order 1/ J (see Fig. 3). For very small 
J, this negative portion of C{t) almost completely can- 
cels the short-time positive portion, resulting in an ex- 
tremely small K* (see inset in Fig. [3]). This cancelation 
is a hallmark of strong localization and can be observed, 
e.g. in an Anderson insulator where it is nearly com- 
plete {k* exponentially with time). While the very 
short time behavior at small J is easily reproduced an- 
alytically by ignoring dynamical spin-spin correlations, 
the behavior out to times of order 1/J is representative 
of correlated motion of few spins (likely pairs) . Although 
likely non-integrable, this motion is nevertheless mostly 
quasipcriodic — we recorded indications of this in local 
spin-spin correlation functions (not shown here). 

Finally, there is apparently a power-law long-time tail 
with a negative amplitude: C{t) ^ — t^^^^, with an ex- 
ponent that we find is approximately x = 0.25 over an 
intermediate range of 0.2 ^ J ^ 0.4 (and more generally. 




FIG. 2: (color online) Short time behavior of C{t) for J = 0.32 
(red, noticeably different trace) and J — 0.08, 0.12, 0.16 (these 
are almost identical data in this plot). Note rescaling of the 
vertical axis by J^. 

0.005 r 




FIG. 3: (color online) Current autocorrelations on medium 
time scales ~ 1/J for J = 0.32, 0.16, 0.12, 0.08, from top (red) 
to bottom (green) trace at tJ = 2). Note the rescaling of both 
the vertical and time axes. The inset shows near cancellation 
between short and medium times. 

perhaps). To observe this with the least amount of effort 
it is best to average C{t) at long times over a neighbor- 
hood of t (see Fig. S]) or to measure KL^t) and compute its 
"exponential derivative", ri{t) = HL{t) ~ KL{2t), at a se- 
quence of points i„ = 2", n = 1, 2, 3 . . . (see Fig. [SJ. The 
apparent value x = 0.25 of this exponent is something 
that we do not understand yet theoretically. However, 
we find that it does provide a good fit to the data over a 
wide dynamic range, providing some support for our use 
of it to extrapolate to infinite time and thus estimate the 
DC thermal conductivity, as discussed below. 



B. DC conductivity: extrapolations and fits 

Our extrapolations of the DC conductivity will be 
based entirely on the long time behavior of kl (t) evalu- 
ated at a set of times t„ = 2" with integer n and for large 
enough L to eliminate finite-size effects (so we drop the 
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FIG. 4: (color online) Top panel: long-time tail in the current 
autocorrelation function for J = 0.20,0.24,0.28,0.32,0.40 
shown bottom to top in red, orange, yellow, green and blue, 
respectively. Bottom panel; to estimate the exponent we mul- 
tiply the data by t^''* (and also display lines with slope ±0.05). 
Although these data do not exclude an exponent that varies 
with J, we interpret these results as supportive of a single 
exponent x « 0.25 at asymptotically long times but with a 
more pronounced short-time transient at smaller J. 
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FIG. 5: (color online) Long time tails as seen from ri{t) 
for J = 0.20,0.24,0.28,0.32,0.36,0.40,0.48 (bottom to top). 
Black line is a guide to the eye with slope —1/4. Note that 
the short-time transients are stronger here, as compared to 
the auto-correlation data in Fig. 4. 



subscript L). We start by describing the procedure used 
to arrive at the numerical estimates of the DC conduc- 




FIG. 6: Variation of K{t) for J = 

0.64, 0.56, 0.48, 0.40, 0.36, 0.32, 0.28, 0.24, 0.20, 0.18, 0.16, 0.14, 0.12 
plotted vs. s = (1 -|- Jt)""'^^. Lines are merely guides to the 
eye, and statistical errors are too small to be seen on most 
of these points. This figure is used for obtaining J > 0.36 
entries in Table [T] 

tivity, then turn to the subject of uncertainties. 

A typical instantaneous value of the energy current is 
set by the strength of the exchange, J. As a consequence 
K{t) ~ for small J at short and intermediate times {t 
of order 1/ J or less) . Given the time-dependence at in- 
termediate and long times, as discussed above, we adopt 
the variable s = (1 -|- Jt)~^^^^ as a convenient "scaling" 
of time for displaying our results. These rescalings "col- 
lapse" the observed values of K{t) for short to intermedi- 
ate times across the entire range of J studied, as shown 
in Fig. ini 

The extrapolated values of the DC conductivity de- 
crease strongly as J is reduced. Extrapolation of K{t) to 
s = and thus DC is fairly unambiguous for J > 0.32, as 
can be seen in Fig. [S] To display the long-time results at 
smaller J, in Fig. [7] we instead show k/ J^^. Here one can 
see that as we go to smaller J the extrapolation to the 
DC limit (s = 0) becomes more and more of "a reach" as 
J is reduced. The outcomes of these extrapolations and 
rough estimates of the uncertainties are summarized in 
Table H 

There are several sources of uncertainty in the esti- 
mates of the DC thermal conductivity k reported in Ta- 
ble I. These can be separated into those originating with 
the measured values of kl(^) and those due to the ex- 
trapolation to DC. 

The statistical uncertainties in the measured values of 
KL{t) were estimated (and shown in the figures) from 
sample-to-sample fluctuations which we find follow gaus- 
sian statistics to a good approximation for these long 
(large L) samples. We did look for a possible systematic 
source of error originating with roundoff and its amplifi- 
cation by chaos (see Appendix B) and found it not to be 
relevant for the values of J and t studied. 

The uncertainties in our estimates of the DC k from the 
extrapolation procedure begin with the assumed value of 
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FIG. 7: (color online) Same data as in Fig. 6, but now scaled 
and displayed in a way that allows one to see the extrapo- 
lations to s = {t —> oo) for small J. Note the diflterent 
rescaling schemes used in preceeding and current plots to fo- 
cus on "collapse" of short time data (previous plot) vs. long 
time extrapolations (present plot). As before black lines are 
drawn through the data for guiding the eyes. Colored lines 
are results of polynomial fits, as explained in the text. This 
figure is used for obtaining entries in Table |I] for J < 0.32. 
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500 


27 
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TABLE I: Extrapolated estimates of the DC conductivity k, 
estimated uncertainties, length L of samples, and the number 
of time steps T of the runs. 



the long-time powerlaw, x = 0.25. Clearly, using a dif- 
ferent exponent will change the extrapolated DC values 
of K somewhat. This uncertainty increases with decreas- 
ing J as the ratio of the ^^(i) at the last time point 
to the extrapolated value increases. At our smallest J 
values, the curvature in our k vs. s plots due to the 
crossover to the earlier-time insulating- like ^ 1/t 
dependence becomes more apparent and further compli- 
cates the extrapolation. Although we have experimented 
some with different schemes for extrapolating to DC, in- 
cluding different choices of exponent x, in the end the 
following procedure appeared to capture the overall scale 
of the diffusion constant, and with a generous estimate of 
the uncertainty: i) we start by removing early data with 
s > 0.5 to focus strictly on the long-time behavior; ii) this 
long time dependence is further truncated by removing 
5 latest points and then fitted to a polynomial X^o'^n'S" 
to better capture the curvature apparent in the data - 
these fits are shown in Fig. [7] and qq are the DC values 
reported in Table HI iii) the uncertainty is estimated as 
the greater of statistical error in the last point (which 
is negligible for most of our data) and the difference be- 
tween Oo and a simple linear extrapolation performed on 
latest five data points not included in (ii). 

Overall, we deem the values presented in Table [T] as 
"safe" since all extrapolated k's differ by at most a factor 
2 from k's actually measured, in other words our extrapo- 
lations are reasonably conservative (with the exception of 
two smallest J's where the extrapolations yields stronger 
reductions). 



IV. SUMMARY, FURTHER EXPLORATIONS 
AND OUTLOOK 

In summary, we considered a rather generic model 
of classical Hamiltonian many-body dynamics with 
quenched disorder, and explored the systematic variation 
in the thermal diffusivity between conducting and insu- 
lating states. We found a rapid variation of the diffusion 
constant and presented quantitative estimates of the lat- 
ter across more than 5 orders of magnitude of change. 
The origin of this behavior may be traced to spatial lo- 
calization of classical few-body chaos - we plan to present 
further results along these lines separately. Qualitatively, 
such a scenario is rather plausible at very low J, where 
most spins are spectrally decoupled due to disorder and 
essentially just undergo independent Larmor precessions. 
As long as J is nonzero, however, there will always be a 
fraction of spins in resonance with some of their immedi- 
ate neighbors. These clusters are then deterministically 
chaotic and thus generate broad-band noise, which allows 
them to exchange energy with all other nearby spins. Im- 
portantly, in the entire parameter range studied this het- 
erogeneous regime eventually gives way at long time to 
a more homogeneous conducting state in the DC limit. 
Thus, we suspect that internally generated but localized 
noise always causes nonzero DC thermal transport even 



8 



in the strongly disordered regime, as long as the spin-spin 
interaction J is nonzero. 

Additionally, we also discovered and characterized an 
apparent, novel finite-time (frequency) correction to dif- 
fusion, with the diffusivity varying as D{ijj) k, D{0) + 
a\uj\^ with X = 0.25. Previous theoretical work on cor- 
rections to diffusion due to quenched disorder [3] have in- 
stead found a correction with exponent x = 1/2, which 
is quite inconsistent with our numerical results. This 
powerlaw behavior is apparently not due to the localiza- 
tion of chaos discussed above, as it persists well into the 
strongly conducting regime (larger J) and also exists in 
models without a strong disorder limit at all (e.g. with 
random fields of equal magnitude but random direction; 
data not shown) . So far we have not found a theoretical 
understanding of these interesting corrections to simple 
thermal diffusion. 
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APPENDIX A: FINITE SIZE EFFECTS 

We have checked that all of the extrapolations above 
are free from finite size effects (by comparing against sim- 
ulations on smaller, and in some cases, larger samples). 
Nevertheless, it is interesting to consider the expected 
hydrodynamic size effects somewhat quantitatively, via 
Eq. [151 To illustrate this we display in Fig. [8] some re- 
sults on shorter systems for J = 0.16,0.32,0.40 that do 
show size effects: due to the periodic boundary condi- 
tions, the conductivity in smaller rings saturates in the 
DC limit at a value corresponding to the AC value at 
a "frequency" 27r/r corresponding to = CdD{2tt/L)'^, 
with D = 3k. Our results are qualitatively consistent 
this with Cd — 10 or slightly larger. It perhaps remark- 
able that despite orders of magnitude of variation in the 
diffusion constant in going from J = 0.4 to J = 0.16 the 
crossover from bulk to finite system behavior is charac- 
terized by roughly the same constant = 10. 

APPENDIX B: CHAOS AMPLIFICATION OF 
ROUND-OFF ERRORS 

No numerical study of a nonlinear classical dynamical 
system is complete without some understanding of the in- 
terplay of discretization and round-off errors and chaos. 




FIG. 8: (color online) Finite size effects for J = 
0.16,0.32,0.40: 100 vs. 20 spins for J=0.40, 10 and 20 vs. 
100 spins for J=0.32, 20 vs. 100 spins for J=0.16. Red color 
is used to indicate the data influenced by finite size effects 
according to Eq. [TS]with Cd ~ 10. Inset: J = 0.16 data 
replotted. 



We are studying a Hamiltonian system that conserves to- 
tal energy, so the chaos is only within manifolds of con- 
stant total energy in configuration space. Thus although 
round-off errors introduce tiny violations of energy con- 
servation, these changes in the total energy are not sub- 
sequently amplified by the system's chaos; we have nu- 
merically checked that this is indeed the case. As a result 
of this precise energy conservation the energy transport 
computation remains well-defined. The simulation is far 
less stable within an equal-energy manifold, where nearby 
trajectories diverge exponentially due to chaos. In par- 
ticular, this means that the component of any round- 
off error that is parallel to the equal-energy manifold is 
exponentially amplified by the chaos. At large J this 
happens rather quickly, while for small J the chaos is 
weaker and longer individual trajectories can be retraced 
back to their respective initial conditions. However, at 
small J very long runs are necessary to extrapolate to 
the DC thermal conductivity: in the end all of our ex- 
trapolations are done in the regime where all individual 
many-body trajectories are strongly perturbed by chaos- 
amplified round-off errors. 

Ultimately, however, we are only concerned with the 
stability of the current autocorrelations that enter in the 
Kubo formula for k. Although the precise trajectories 
may diverge due to chaos-amplified round-off errors, this 
need not have a strong effect on C{t). To study this issue 
quantitatively we simulated roundoff noise of different 
strength in our computations. Specifically, we add extra 
random noise to the computation without altering the 
total energy by multiplying the angle each spin precesses 
in each time step by a factor of 1 -I- ?7i(i), where the rii{t) 
are independent random numbers uniformly distributed 
between P and —P [P = noise strength). 

In 400 rings of 500 spins coupled with J = 0.14 
we simulated the same initial condition with different 
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0.07 r 




FIG. 9: Roundoff effects at J = 0.14: low frequency, long- 
time conductivity is larger for larger values of P but is essen- 
tially indistinguishable between P = and P = 0.01. Inset: 
data with P = 0.1, 0.01,0. 



r]i{t), and with different values of simulated noise P = 
10°, 10~^, 10"^, - these results are presented in Fig. [5] 
below. As expected, the long-time insulating behavior 
is weakened by the presence of noise. Quantitatively, 
however, we observe little or no difference between re- 
sults obtained in the presence of simulated noise with 
P = 10~^ vs. ones obtained for intrinsic noise (which at 
double precision corresponds to Pmtr ^ 10^^^). Clearly, 
this statement heavily depends on the duration of the 
simulation, value of J, etc. Judging from Fig. [5] round- 
off errors are not a serious source of uncertainty in our 
results in main text. Interestingly, it is also possible for 
strong noise to suppress k, as indeed happens at shorter 
times, which can be traced here to a sort of "dephasing" 
of sharp response of quasiperiodic localized states. 
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